home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / eigenlib / hqr.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-12-15  |  7.3 KB  |  223 lines
  1. hqr(nm,n,low,igh,h,wr,wi)
  2.  
  3. int n,nm,*igh,*low;
  4. double **h,*wr,*wi;
  5. {
  6. int i,j,k,l,m,en,ll,mm,na,itn,its,mp2,enm2,ierr;
  7. double p,q,r,s,t,w,x,y,zz,norm,tst1,tst2;
  8. double fabs(),sqrt(),copysign();
  9. char notlas;
  10. /*
  11.       this subroutine is a translation of the algol procedure hqr,
  12.       num. math. 14, 219-231(1970) by martin, peters, and wilkinson.
  13.       handbook for auto. comp., vol.ii-linear algebra, 359-371(1971).
  14.  
  15.       this subroutine finds the eigenvalues of a real
  16.       upper hessenberg matrix by the qr method.
  17.  
  18.       on input
  19.  
  20.          nm must be set to the row dimension of two-dimensional
  21.            array parameters as declared in the calling program
  22.            dimension statement.
  23.  
  24.          n is the order of the matrix.
  25.  
  26.          low and igh are integers determined by the balancing
  27.            subroutine  balanc.  if  balanc  has not been used,
  28.            set low=1, igh=n.
  29.  
  30.          h contains the upper hessenberg matrix.  information about
  31.            the transformations used in the reduction to hessenberg
  32.            form by  elmhes  or  orthes, if performed, is stored
  33.            in the remaining triangle under the hessenberg matrix.
  34.  
  35.       on output
  36.  
  37.          h has been destroyed.  therefore, it must be saved
  38.            before calling  hqr  if subsequent calculation and
  39.            back transformation of eigenvectors is to be performed.
  40.  
  41.          wr and wi contain the real and imaginary parts,
  42.            respectively, of the eigenvalues.  the eigenvalues
  43.            are unordered except that complex conjugate pairs
  44.            of values appear consecutively with the eigenvalue
  45.            having the positive imaginary part first.  if an
  46.            error exit is made, the eigenvalues should be correct
  47.            for indices ierr+1,...,n.
  48.  
  49.          ierr is set to
  50.            zero       for normal return,
  51.            j          if the limit of 30*n iterations is exhausted
  52.                       while the j-th eigenvalue is being sought.
  53.  
  54.  
  55.       this routine is a C-translation of the FORTRAN 77 source code
  56.       written by the mathematics and computer science division,
  57.       argonne national laboratory
  58.       last change :   september 1989.
  59.  
  60.       mark myers
  61.       Center for Applied Mathematics 
  62.       Cornell University    (607) 255-4195
  63.       --------------------------------------------------------- */ 
  64.  
  65.         ierr = 0;
  66.         norm = 0.e0;
  67.         k = 1;
  68.         for (i=1;i<=n;i++)                   /* store roots isolated by balanc */
  69.           {for (j=k;j<=n;j++)                /* and compute matrix norm.       */
  70.              norm = norm + fabs(h[i][j]);
  71.            k = i;
  72.            if (i < *low | i > *igh)
  73.              {wr[i] = h[i][i];
  74.               wi[i] = 0.e0;}  }
  75.         en = *igh;
  76.         t = 0.e0;
  77.         itn = 30*n;
  78. cycle1:  ;
  79.         if (en < *low)
  80.           return;
  81.         its = 0;
  82.         na = en - 1;
  83.         enm2 = na - 1;
  84. cycle2: ;
  85.         for (ll = *low;ll <= en;ll++)    /* look for single small sub-diagonal element  */
  86.           {l = en + *low - ll;      /* for l=en step -1 until low do               */
  87.            if (l == *low) 
  88.              break;
  89.            s = fabs(h[l-1][l-1]) + fabs(h[l][l]);
  90.            if (s == 0.e0)
  91.          s = norm;
  92.            tst1 = s;
  93.            tst2 = tst1 + fabs(h[l][l-1]);
  94.            if (tst2 == tst1) 
  95.         break;}
  96.  
  97.         x = h[en][en];                 /* form shift */
  98.         if (l == en)
  99.           {wr[en] = x + t;
  100.            wi[en] = 0.e0;
  101.            en = na;
  102.            goto cycle1;}
  103.         y = h[na][na];
  104.         w = h[en][na] * h[na][en];
  105.         if (l == na) 
  106.           {p = (y - x) / 2.e0;
  107.            q = p * p + w;
  108.            zz = sqrt(fabs(q));
  109.            x = x + t;
  110.        if(q<0.e0)
  111.              {wr[na] = x + p;           /* complex pair */
  112.               wr[en] = x + p;
  113.               wi[na] = zz;
  114.               wi[en] = -zz;
  115.               en = enm2;
  116.               goto cycle1;}
  117.            zz = p + copysign(zz,p);        /* real pair */
  118.            wr[na] = x + zz;
  119.            wr[en] = wr[na];
  120.            if (zz != 0.e0)
  121.          wr[en] = x - w / zz;
  122.            wi[na] = 0.e0;
  123.            wi[en] = 0.e0;
  124.        en = enm2;
  125.        goto cycle1;}
  126.         if (itn == 0) 
  127.           return(en);
  128.         if (its == 10 | its == 20)
  129.          {t = t + x;                  /* form exceptional shift */
  130.           for (i = *low;i<=en;i++)
  131.             h[i][i] = h[i][i] - x;
  132.           s = fabs(h[en][na]) + fabs(h[na][enm2]);
  133.           x = 0.75e0 * s;
  134.           y = x;
  135.           w = -0.4375e0 * s * s;}
  136.         its = its + 1;
  137.         itn = itn - 1;
  138.         for (mm=l;mm<=enm2;mm++)       /* look for two consecutive small */ 
  139.           {m = enm2 + l - mm;         /* sub-diagonal elements.         */
  140.            zz = h[m][m];              /* for m=en-2 step -1 until l...  */
  141.            r = x - zz;
  142.            s = y - zz;
  143.            p = (r * s - w) / h[m+1][m] + h[m][m+1];
  144.            q = h[m+1][m+1] - zz - r - s;
  145.            r = h[m+2][m+1];
  146.            s = fabs(p) + fabs(q) + fabs(r);
  147.            p = p / s;
  148.            q = q / s;
  149.            r = r / s;
  150.            if (m == l)
  151.          break;
  152.            tst1 = fabs(p)*(fabs(h[m-1][m-1]) + fabs(zz) + fabs(h[m+1][m+1]));
  153.            tst2 = tst1 + fabs(h[m][m-1])*(fabs(q) + fabs(r));
  154.            if (tst2 == tst1)
  155.          break;}
  156.         mp2 = m + 2;
  157.  
  158.         for (i=mp2;i<=en;i++)
  159.           {h[i][i-2] = 0.e0;
  160.            if (i == mp2)
  161.              continue;
  162.            h[i][i-3] = 0.e0;}
  163.         for (k=m;k<=na;k++)        /* double qr step involving rows l to en & */
  164.           {if (k==na)             /* columns m to en                         */
  165.              notlas = 'f';
  166.            else
  167.            notlas = 't';
  168.            if (k != m) 
  169.              {p = h[k][k-1];
  170.               q = h[k+1][k-1];
  171.               r = 0.e0;
  172.               if (notlas == 't')
  173.               r = h[k+2][k-1];
  174.               x = fabs(p) + fabs(q) + fabs(r);
  175.               if (x == 0.e0) 
  176.               continue;    
  177.               p = p / x;
  178.               q = q / x;
  179.               r = r / x; }
  180.            s = copysign(sqrt(p*p+q*q+r*r),p);
  181.            if (k != m) 
  182.              h[k][k-1] = -s * x;
  183.        else 
  184.              if (l != m)
  185.            h[k][k-1] = -h[k][k-1];
  186.            p = p + s;
  187.            x = p / s;
  188.            y = q / s;
  189.            zz = r / s;
  190.            q = q / p;
  191.            r = r / p;
  192.            if (notlas != 't')
  193.              {for (j=k;j<=n;j++)                      /* row modification */
  194.                {p = h[k][j] + q * h[k+1][j];
  195.                 h[k][j] = h[k][j] - p * x;
  196.                 h[k+1][j] = h[k+1][j] - p * y;}
  197.    
  198.             j = en;
  199.             if(en>k+3)
  200.               j = k+3;
  201.               for (i=1;i<=j;i++)                     /* column modification */
  202.                 {p = x * h[i][k] + y * h[i][k+1];
  203.                  h[i][k] = h[i][k] - p;
  204.                  h[i][k+1] = h[i][k+1] - p * q;}}
  205.            else      
  206.              {for (j=k;j<=n;j++)                     /* row modification */
  207.                {p = h[k][j] + q * h[k+1][j] + r * h[k+2][j];
  208.                 h[k][j] = h[k][j] - p * x;
  209.                 h[k+1][j] = h[k+1][j] - p * y;
  210.                 h[k+2][j] = h[k+2][j] - p * zz;}
  211.     
  212.           j = en;
  213.           if (k+3<en)
  214.             j = k+3;
  215.               for (i=1;i<=j;i++)                    /* column modification */
  216.                 {p = x * h[i][k] + y * h[i][k+1] + zz * h[i][k+2];
  217.                  h[i][k] = h[i][k] - p;
  218.                  h[i][k+1] = h[i][k+1] - p * q;
  219.                  h[i][k+2] = h[i][k+2] - p * r;}}
  220.     }
  221.       goto cycle2;
  222. }
  223.